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ABSTRACT 

Using a high resolution dark matter (DM) simulation of the Local Group, conducted 
within the framework of the Constrained Local UnivcrsE Simulation (CLUES) project, 
we investigate the nature of how satellites of the Milky Way (MW) and M31 are 
accreted. Satellites of these two galaxies are accreted anisotropically onto the main 
halos, entering the virial radius of their hosts, from specific "spots" with respect to the 
large scale structure. Furthermore, the material which is tidally stripped from these 
accreted satellites is also, at z = 0, distributed anisotropically and is characterized by 
an ellipsoidal sub-volume embedded in the halo. The angular pattern created by the 
locus of satellite infall points and the projected z = stripped dark matter stripped is 
investigated within a coordinate system determined by the location of the Local Group 
companion and the simulated Virgo cluster across concentric shells ranging from 0.1 
to 5r v i r - Remarkably, the principal axis of the ellipsoidal sub- volume shows a coherent 
alignment extending from well within the halo to a few r v j r . A spherical harmonics 
transform applied to the angular distributions confirms the visual impression: namely, 
the angular distributions of both the satellites entry points and stripped DM for 
both halos is dominated by the / = 2 quadrupolc term, whose major principal axis 
is approximately aligned across the shells considered. It follows that the outer (r > 
0.5r v ; r ) structure of the main halos of the Local Group composed of stripped material is 
closely related to the cosmic web, within which it is embedded. Given the very plausible 
hypothesis that an important fraction of the stellar halo of the Milky Way has been 
accreted from satellite galaxies, the present results can be directly applied to the stellar 
halo of the MW and M31. We predict that the remnants of tidally stripped satellites 
should be embedded in streams of material composed of dark matter and stars. The 
present results can therefore shed light on the existence of satellites embedded within 
larger streams of matter, such as the Segue 2 satellite. 



INTRODUCTION 



In the ACDM model, Galaxy halos like that encompass- 
ing our own Milky Way, are fo rmed through the accre- 
tion of smaller substructures (e.g. IWhite fe R ccs 1978]). Nu- 
merical simulations of the formation of galaxy sized halos 
have consistently shown that small, self bound, over-dense 
clumps of matter can survive the violent physical processes 
associated with tidal stripping and dynamical friction and 
end up as distinct structures embedded within larger halos 
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Once accreted, some of these so-called substructures are nat- 
urally identified as the sites of satellite galaxies and their 
properties can be thus be compared to observations of satel- 
lites in our own halo. 

The satellites of the Milky Way's halo posit a number of 
interesting questions. That they are all aligned on a great cir- 
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cle in the sky has l ong been curi ous (IKunkel fc Demersll 19761 ; 
iLvnden-Belll 19821 ). Yet recentlv lKroupa et all (|2005l ) argued 



that the anisotropic distribution of the Milky Way's 11 "clas- 
sical" satellites is inconsistent with the hypothesis that lu- 
minous satellites are a random sub-sample of dark matter 
substructures. Furthermore the newly discov ered Sloan Dig- 
ital S ky Survey (SDSS) satellites (e.g. see iKoposov et al.l 
120081 . and references therein) appear to be loca ted directly 
on this great disc-of-s atellites (iMetz et al.ll2009j). 

Many studies dLibeskind et all 120051 . 12001 |2009| ; 
IZentner et alj 120051 ; iKang et al. 20051 ) have employed nu- 
merical si mulations in a b id to understand this peculiar 
geometry. iLibeskind et al.l (|2005l ) found that the flattened 
anisotropic distribution of satellites is readily reproduced if 
the luminous satellites of the Milky Way are associated with 
substructures that had the largest (pre-infall) progenitors. 
These subhaloes were born in the dense spines of the fil- 
ament that co l lapsed to form the aspherical galactic halo. 
IZentner et alj (|2005l ) showed that the subhalo population 
(albeit in Galactic haloes simulated to sit within filaments) 
as a whole is not isotropically distributed but rather tend to 
occupy sites close to the long axis of the parent halo. They 
thus argued that an isotropic subhalo population is not the 
correct null-hypothesis for the ACDM model 

While the planarity of Milky Way satellites is no longer 
deemed a threat to the standard model, its origin has evaded 
a definitive understanding. Most authors have indeed man- 
aged to reproduce the observed anisotropy and have only 
hinted at what could be its cause. One of the most inter- 
esting ideas is that th e satellites o f the M ilky Way entered 
the halo as a group. iLi fc Helmil (|2008l ) analyzed a suite 
of high resolution dark matter only simulations populated 
with galaxies semi-analytically and found that around one 
third of all satellites were accreted as part of larger groups 
of satellites. This result is promising since it implies a flat- 
tened z = distribution of satell ites. However in cons t rained 
simulation of the Local Group, iKlimentowski et alj (|2009T ) 
found that although satellites form in groups, they tend to 
be loosely bound, and the group identity is short lived. 

A complimentary explanation of the flattened distribu- 
tion of satellites is that subhaloes are not uniformly accreted 
from all directions in the sky. Instead, substructures stream 
across intergalactic space through filaments that feed the 
growth of a Galaxy-sized halo w hich has been exp licitly 
shown for Cluster-sized haloes by iKnebe" et all (|2004l ). Yet 
often the effect of the large scale environment is ignored in 
studies of the spatial distribution of satellite galaxies. 
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initial condition^ that aim at reproducing the particulars of 
the observed Local Group in order to study if substructures 
enter the halo with a uniform distribution with respect to the 



1 See jhttp: / /wwwT CLUES-project.org for a summary of the Cc 
strained Local UnivcrsE Simulation Project. 



large scale structure or whether they come from preferred 
directions. We examine the effect preferred infall has on the 
stripped material that is brought in by accreted subhalos. 

This paper is organized as follows: in Sec.[5]we introduce 
our simulation, the subhalo finder that was used, and the 
large scale geometry. In Sec. [3] we present our results and we 
conclude in Sec. U 



2 METHODS 

In this section we describe in brief the simulations used as 
well as the halo and subhalo finding algorithm employed to 
identify satellites. 



2.1 Constrained Simulations of the Local Group 

The simulations used in this work are embedded in the 
Constrained Local UniversE Simulation (CLUES) project 
and have been already studied i n a number of recent pa- 
pers (e.g Libeskind et alj 120091 ; IKlimentowski et all 120091 ; 
IKnebe et alj 20091). We refer th e reader to those papers (in 
particular Libeskind et alj|2009l ) for details on how the con- 
straints were generated and how the simulations were run: 
we highlight just the salient points here for clarity. 

We choose to run our simulations using standard ACDM 
initial conditions, that assume a Wilkinson Micr owave 
Anisotropy Probe 3 cosmology (|Spergel et alj 120071 ) , i.e. 
ilm = 0.24, Q b = 0.042, D.a = 0.76 and h = 0.73. We use 
a normalization of <rg = 0.73 and an n = 0.95 slo pe of the 
powe r spectrum. We used the MPI code gadget2 (|Springell 
l2005h to simulate the evolution oft a cosmological box with 
side length of Lbox = 64h, _1 Mpc. 

Instead of seeding our initial conditions as a just a ran- 
dom cube of space, the initial conditions of our volume are 
constrained to reproduce, at 2 = a number of objects that 
compose the local environment (see lHoffman fc Ribak|[l99ll . 
for details on how constraints initial conditions are gener- 
ated), including a "virgo" cluster, a "coma" cluster and a 
"Local Group" . Our method allows us to properly constrain 
the large scales (i.e. those still linear by 2 = 0) but we do 
not constrain the local group itself. In order to obtain a lo- 
cal group in the correct environment, three low resolutions 
constrained simulations are run with varying random seeds. 
Each 2 = low resolution simulations is then examined 
and if an object that resembles the local group is found (by 
construction this will be in the correct place), these initial 
conditions are selected for high resolution re-simulation. 

We resimulate just the region of interest around the 
local group. We center a sphere of radius 2 /i _1 Mpc around 
the local group and populate with ~ 5.2 x 10 7 low mass, high 
resolution particles. Within our local group we are thus able 
to achieve a particle mass of just M dm = 2.54 x 10 5 /i _1 M . 

Our constraints reproduce a cosmography which closely 
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Property 


Simulated LG 


Observed LG 


Reference 


M MW 


6.57 x 10 n M Q 


1O 12 M 


[1,2,3] 


MM31 


8.17 x 10 11 AI e 


8.2 x 10 n M Q 


[4] 


Myirgo 


1.08 x 1O 15 A1 


1.2 x 1O 15 M 


[5] 


^MW-M31 


0.909 Mpc 


0.77 Mpc 


[6,7] 


^MW-Virgo 


15.1 Mpc 


16.5 Mpc 


[8] 


Vm31 


142 kms" 1 


195 kms- 1 


[9] 



Table 1. The 2 = properties of the simulated and observed 
Local Group. From the top row own, we shoe the following prop- 
erties: the mass of the MW's halo (Mjyiw), the mass of M31's halo 
(Mm3i), the mass of the Virgo cluster (Mvirgo), the distance from 
the MW to M31 (cZmw-M3i), the distance from the MW to Virgo 
(gJmw— Virgo), and the relative velocity of M 31(Vm3i ). The refer - 
en ces are as follows: [ljIXue et al.1 (120081') ; [2llKlypin et alj J2002I); 
[31 ISmith et al.l [1200711; [41 ISeigar et alj [|2008h;[5l IFouqu£ et all 
112001 ); [6] Karachcntscv_et_aL I j2004l); [71 IVilardell et al.l ll2010t) ; 



[8l lMei et al.1 d200^ [91 Ivan der Marel fc Guhathakurtal (I2008T) 



resembles the observed Local Group. In Table. [T]we compare 
properties of the simulated local group with observations of 
the real on^l- Although our results do not match the ob- 
servations perfectly, the cosmography simulated using our 
constraints captures the essence - in terms of mass and dis- 
tances - of the observed Local Group. 



2.2 The halo and subhalo finding algorithm 

In this section, we explain how our halo and subhalo finding 
algorithm works. In order to identify halos and subhaloes 
in our simulation we have run the MPI+OpenMP hybrid 
halo finder AHF (AMIGA halo finder, to be downloaded freely 
fr om http://popia.ft . uam.e s/AMIGA) described in detail 
in lKnollmann fe Knebej (120091) . A HF is an improvement of the 



MHF halo finder ( Gill et al. 20041 ) . which locates local over- 



densities in an adaptively smoothed density field as prospec- 
tive halo centers. The local potential minima are computed 
for each of these density peaks and the gravitationally bound 
particles are determined. Only peaks with at least 20 bound 
particles are considered as haloes and retained for further 
analysis. In practice for this work, we only consider subhalos 
with more than 100 particles. We would like to stress that 
our halo finding algorithm automatically identifies haloes, 
sub-haloes, sub-subhaloes, etc. For more details on the mode 
of operation and actual fun ctionality we refer the read er to 
the code description paper (|Knollmann fc Knebej[2009l ). 

For each halo, we compute the virial radius r v i r , that is 
the radius r at which the density M(< r)/ (4nr 3 /3) drops 
below AvirPback- Here Pback is the cosmological background 
matter density. The threshold A v i r is computed using the 

2 In a future paper we intend to study in detail the cosmography 
produced by our constrained simulations. 



spherical top-hat collapse model and is a function of both 
cosmological model and time. For the cosmology that we are 
using, Avir = 355 at z = 0. 

Subhaloes are defined as haloes which lie within the 
virial radius of a more massive halo, the so-called host halo. 
As subhaloes are embedded within the density of their re- 
spective host halo, their own density profile usually shows a 
characteristic upturn at a radius rt < rvi r , where r v [ T would 
be their actual (virial) radius if they were found in isolation^ 
We use this "truncation radius" r t as the outer edge of the 
subhalo and hence subhalo properties (i.e. mass, density pro- 
file, velocity dispersion, rotation curve) are calculated using 
the gravitationally bound particles inside the truncation ra- 
dius rt- For a host halo we calculate properties using the 
virial radius r v i r . 

We build merger trees by cross-correlating haloes in 
consecutive simulation outputs. For this purpose, we use 
a tool that comes with the AHF package and is called 
Merger/Tree. As the name suggests, it serves the purpose 
of identifying corresponding objects in the same simulation 
at different redshifts. We follow each halo (either host or 
subhalo) identified at redshift z — backwards in time, 
identifying as the main progenitor (at the previous redshift) 
the halo that both shares the most particles with the present 
halo and is closest in mass. The latter criterion is important 
for subhaloes given that all their particles are also typically 
bound to the host halo, which is typically orders of magni- 
tude more massive. Given the capabilities of our halo finder 
AHF and the appropriate construction of a merger tree, sub- 
haloes will be followed correctly along their orbits within the 
environment of their respective host until the point where 
they either are tidally destroyed or directly merge with the 
host. 

In practice, we use a lower particle limit of 100 par- 
ticles when considering subhalos. In order to ensure that 
our results are not resolution dependent, we have also run 
our analysis using only subhalos with > 500 particles. The 
results are fully consistent with those presented here, albeit 
with the disadvantage of considering far fewer subhalos. Our 
particle limit constrains us to only study the haloes history 
after z ~ 2, which corresponds to more than 10 Gyrs ago. 
This is because before this redshift, the majority of accreted 
objects fall below our particle limit. Lowering the particle 
limit used to consider subhalos would have the added advan- 
tage of allowing us to probe earlier times, but would return 
subhalos that were not fully resolved. Thus in the rest of 
the text when we refer to the "entire history of the halo" , we 
mean since z ~ 2. Furthermore, by z = 2 the MW and M31 

3 Please note that the actual density profile of subhaloes after the 
removal of the host's background d rops faster than for isolated 
haloes (e.g. LCazantzidis et al. 2 004T) ; only when measured within 
the background still present will we find the characteristic upturn 
used here to define the truncation radius rt. 
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halos have assembled just 7% and 5% of their final mass re- 
spectively, indicating that the z < 2 time period corresponds 
to the main formation epoch of our halos and that the period 
z > 2 is negligible from a mass accretion perspective. 



3 RESULTS 

In this section we present our results regarding the preferen- 
tial infall direction of subhalos. We begin by describing the 
coordinate system we use in order to fix the axes from which 
we measure the infalling direction. This coordinate system 
serves as the basis for determining the preferred direction of 
accretion. 



3.1 The Milky Way - Andromeda - Virgo 
coordinate system 

Since our constraints accurately simulate the local environ- 
ment, we define a coordinate system which is fixed with 
respect to the large scale structure. We show this coordi- 
nates system in Fig. [1] From our Merger tree, we know the 
location of the two halos we choose to call t he "Milky Way" 
(MW ) and "Andromeda" (M31 ) (e.g. see iLibeskind et all 
120091 ; iKlimentowski et al.l 120091 ) at each snapshot. Addi- 
tionally, from our constraints, we know the location of the 
"virgo" cluster. We thus define at each snap shot a coor- 
dinate system centered on each of these galaxies. We call 
7?mw-m31, the position vector connecting MW to M31 (and 
vice versa) the 'V-axis. Since the two galaxies move with 
roughly the same velocity in roughly the same direction, this 
vector changes little over time. We then examine -Riviw-vh-go 
and i?M3i -Virgo, the position vector of the Virgo cluster 
(some 10 Mpc away) centered on the MW and M31. We 
dot i?MW-virgo and i?M3i-vir g o with i?MW-M3i in order to 
obtain the component of the i?MW-vir go that is perpendicu- 
lar to i?MW-M3i , calling this vector our "?/"-axis. In practice, 
this vector is (co-incidentally) nearly 90 degrees away from 
i?MW-M3i and so requires little adjustment. Once we have 
an x and y axis, we compute a z axis by crossing x with y 
(i.e. z — x x y). 

The coordinate system, shown in Fig. [1] provides us 
with three axes that at each time-step are recalculated as the 
systems evolves. In practice, the coordinate system changes 
little, as both M31 and the MW move with similar velocities. 
We may use these axes to note the location on the sky each 
and every time a subhalo crosses the virial radius or factor 
thereof. 



3.2 Prefered directions - a qualitative description 

In this section we present a qualitative description of our re- 
sults regarding the preferred directions we have found. We 



describe the infall of satellites, the unique z — distribu- 
tion of dark matter stripped from these satellites and the 
self similar distribution across a large radial extent. In prac- 
tice in order to do identify the point of satellite infall, we 
examine two adjacent snapshots. We then ask the question: 
of those subhalos within some radius at the later snapshot, 
which had progenitors that were outside of this radius one 
time step earlier? Once all the accreted subhalos have been 
determined, we identify the point of entry. 



3.2.1 The preferential infall of satellite galaxies 

In this section we focus on the locations on the sky from 
which subhalos are accreted. We defer to section [331 a more 
rigorous, quantitative analysis. 

In the top panel of Fig. [2] we examine the anisotropic 
accretion of satellites across concentric shells centered on 
the MW. In the interest of space and focus for all the Aitoff 
plots in this paper, we omit the same plots for M31 but em- 
phasize that the results are qualitatively and quantitatively 
very similar. In an Aitoff projection of the sky, we show the 
location of each satellite's infall point integrated over the en- 
tire history of the halo, rotated into our large scale structure 
coordinate system. In practice, we start by first defining four 
concentric shells at r c — 0.2, 1, 3 and 5 r v i r - We then ex- 
amine two chronologically sequential snapshots and identify 
all the (sub)halos which were accreted across each spherical 
surface, plotting them in the upper panels of Fig. [2] 

We note that neither the MW (shown here) nor M31 
(not shown) display a uniform distribution subhalo accre- 
tion points. Satellite infall points are clearly clustered, rather 
than being uniformly distributed on the sky. 

This is our first result: satellites are not accreted by ha- 
los uniformly over the virial radius; rather, they enter the 
halo from preferred directions. In order to determine whether 
this preferred direction is simply dominated by many subha- 
los entering the hal o as a group over a relatively short period 
of cosmic time (e.g. lLi fc HelmfeOOSl) . we have examined the 
locus of entry points of satellites crossing the virial radius 
as a function of epoch. Each epoch shows a very similar en- 
try pattern. Simply put: at each epoch in the halo's history, 
subhalos have a higher tendency to come from one specific 
direction than from many. 

Remarkably, this preferred direction is similar for halos 
crossing shells from r c > 0.2 r v i r to r c = 5 r v i r , reflective 
of the coherent radial infall across this large rangfl There 
are many ways we can identify preferred directions in our 
distributions. In Sec. 13.31 we perform a more sophisticated 
analysis on the angular distribution of infall points. Here 



4 Although the 0.5 r v ; r shell is not shown here, the self similar 
in-fall pattern begins becoming discernible at r Kb 0.5 r v i r . 
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Figure 1. A 2/i _1 Mpc, z = slice through our simulation showing the coordinate system employed. We have rotated the simulation into 
the MW's reference frame. The blue arrows show the coordinate axes used, while the dashed blue line shows the direction to the Virgo 
cluster, some 10 Mpc away. The red arrows (not to scale) indicate the relative velocity of each halo, while the red circles indicate the 
virial radius (to scale). The green lines indicate the preferential infall direction of satellites across r v j r as determined by the eigenvector 
corresponding to the largest eigenvalue of the quadrupole tensor (see sec. I3.3.2)l . In the panel on the left, the black dots correspond to 
a random 0.1% of the dark matter particles. On the right we show a logarithmic density map (linearly contoured) of all the particles in 
the volume. White regions correspond to high projected densities while red and black regions correspond to lower density regions. Note 
that the green lines are roughly aligned in the same direction as the nearly-continuous white region, a filamentary like structure. 



we simply qualify our visual impression by defining a "sin- 
gle" preferred direction, as that infall point with the highest 
number of other infall points within 20° of it. This point is 
indicated in the upper panels of Fig. [2] as the red circle. We 
also plot it in Fig.[T]as the green arrows. We omit this for the 
left most panel (0.2r v j r ), since the distribution is consistent 
with random, the preferred entry point is meaningless. For 
the other three shells at 1, 3 and 5 r v i r the distributions are 
highly clustered in a patch on the sky whose size is roughly 
20° x 20° and at a longitude of between 6 ~ 50° - 70° and 
a latitude between 4> ~ 30° — 50°. We note however, that 
this naive approach tells us little of the strength of the clus- 
tering. Although the filamentary nature of the cosmic web 
can be invoked to explain this affect on large scales, its is 
remarkable that the inflow pattern persists to r ~ 0.5r v ir 
(though not shown here). Note that the preferred directions 
changes little across 5 r v i r . We explore the self-similar nature 
of satellite accretion across many scales in section [3.2.31 

It is not directly clear from the simulations what deter- 
mines this preferred point of entry, although some observa- 
tions can be made. For one, the subhalos tend to originate 
from roughly the same region in the sky for both of the two 



halos: a region in the direction of Virgo. This can be seen in 
Fig Q] as the green arrows. In Fig. [T] right panel, we show a 
logarithmic density map of the volume under consideration, 
in order to highlight the z = DM particle distribution. 
Note that the green arrows can be seen to point in the direc- 
tion of a finger like, filamentary structure. It becomes appar- 
ent from this figure, that as the MW and M31 travel towards 
virgo they scoop up subhalos in their path. These subhalos 
reside in a filamentary like structure that also points in the 
direction of Virgo. As the main halos approach a subhalo, 
the subhalo comes to a "fork in the road" and drifts (per- 
haps randomly) towards one of the parent halos, later to be 
accreted. Identifying the filament in the surrounding matter 
is however, beyond the scope of this paper. 



3.2.2 The z — distribution of dark matter stripped from 
accreted satellites 

In this section we look at the z = distribution of dark 
matter stripped from the preferentially accreted satellites 
described above. 

We wish to see if there is any hint of the preferential 
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shell at 0.2r shell at 1r shell at 3r shell at 5r . 




high density low density 



Figure 2. Top Panels: The point of satellite infall on concentric shells of radius 0.2, 1, 3 and 5r v ; r . Each time a satellite is accreted 
across one of these shells, its entry point is recorded and plotted. The red circle shows the "single" preferred direction, defined as the point 
on the sky with largest number of satellite infall points with 20° of it. Bottom Panels: The projected dark matter density of particles 
brought in by satellites crossing these four radii integrated over the lifetime of the halo. Red regions indicate high density while blue 
regions indicate low density. Regions on the sky where many satellites entered a given shell, correspond to regions of higher projected 
dark matter density. 

onto the surface of a unit sphere and calculate a "density" 
for each particle as the (angular) distance to the 5th nearest 
neighbor. 

We show our results in the bottom panels of Fig. [2] for 
satellites crossing and stripped dark matter particles within 
shells of radius, r c = 0.2, 1, 3, and 5r v i r - A number of salient 
points can be made from this figure. 

Firstly, by comparing the top and bottom panels of 
Fig. there appears to be a strong correspondence between 
the z — angular distribution of stripped dark matter and 
the locus of satellite entry points. Although the correspon- 
dence is not one to one, regions from where many satel- 
lites entered the virial radius coincide with regions of high 
stripped particle density. Likewise, regions where few satel- 
lites entered the halo have correspondingly fewer dark mat- 
ter particles. This brings us to our second main result: dark 
matter that is stripped from subhalos, displays at z = a 
memory of the region from where the parent subhalo was ac- 



infall of satellites, ingrained in the z — distribution of dark 
matter particles. With the exception of the largest subhaloes 
and those accreted a t low z, many accreted subhalos will be 
tidally stripped (e.g. I Johnston 1998; Pc harrubia et al1l2006l ; 
IWarnick et al.l [2009) an d by z = will have lost a large 
amount of their material. We can identify this material by 
noting, at each accretion event, which dark matter particles 
are bound to which subhalo. We then excise from this list, 
the particles which at z = are still bound to clumps. We 
thus obtain a list of dark matter particles that were accreted 
in subhalos and subsequently stripped from them, over the 
entire history of the halo. This list can then be used to locate 
the dark matter particle's z = positions. 

For each halo that is accreted across a given radius (r c ), 
we begin by noting the ID's of its dark matter particles. We 
then locate these particles at z = and consider only those 
within r c . Using the reference frame defined by the large 
scale structure at z — we then project this set of particles 
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creted. We therefore can infer the following picture of sub- 
halo infall: subhalos are accreted from a preferred direction 
over the history of the halo. As these subhalos move through 
the halo, they suffer mass loss effects due to tidal forces. 
Their constituent particles get stripped and the subhalos - 
at some point - cease to exist as recognizable over dense 
clumps. The stripped material however, continues to reside 
in regions centered on where the subhalo was accreted. 

We note that stripped debris material is a sub-set of 
halo particles not bound to substructures. The distribution 
of the latter is mildly aspherical while the debris material 
appears significantly flattened. 

We wish to draw the reader's attention to the self- 
similar patterns across the three shells (where r c > 0.2 r vlI ) 
shown in Fig. [2] A visual inspection immediately reveals 
that the preferred direction on the sky from where satellites 
are accreted changes little from 1 to 5 r v j r . Similarly, the 
2 = dark matter brought in by these satellites also oc- 
cupies similar regions of the sky. Broadly speaking, regions 
of high satellite infall point density correspond to regions of 
high z = stripped dark matter density and these regions are 
constant across an order of magnitude in radial distance. 

We would like to remind the reader that the projected 
angular surface density maps in (the bottom panel of) Fig. 2 
are cumulative, i.e. all the stripped material within the re- 
spective radius is considered. This means that some parti- 
cles may (or may not) contribute to more then one map. 
Consider a subhalo that crosses both 5 and 3 r v i r yet only 
deposits its stripped material interior to 3 r v i r : its stripped 
material will contribute to the density map within both 3r v i r 
and 5r v i r hence possibly enhancing the signal of alignment. 
To gauge the relevance of this "particle stacking" we ex- 
amined surface density maps considering just the stripped 
material that at z = is located within the concentric shells 
probing the regions: 0.2r v i r < r < r v i r , r v - lr < r < 3r v i r , and 
3r v i r < r < 5r v i r - We qualitatively recover the same signal 
that shows the same characteristics as the maps presented 
in the bottom panels of Fig. 2, respectively. 

3.2.3 The self-similar angular distribution of stripped 
material 

In this section we examine in detail the distribution of the 
dark matter stripped from satellites accreted across concen- 
tric shells centered on the galaxy. 

In Fig. [2] we showed that both the locus of satellite infall 
points and the dark matter ellipsoid of stripped particles, 
displays coherence both with each other and across many 
radial scales. We wish to study this ellipsoid in more detail. 
In Fig. we show an Aitoff projection of the z = angu- 
lar position of particles stripped from satellites that crossed 
eight shells at r c = 0.1, 0.2, 0.5, 1, 2, 3, 4, and 5 r v - a . We ro- 
tate each particle distribution into the principal axis defined 
by these particles within r v i r . We find the principal axes of 



the particle distribution by diagonalizing the moment inertia 
tensor: 

Ii.j = ^ m n XiXj (1) 

n 

in order to identify the eigenvectors and eigenvalues (defined 
according to convention as a > b > c). Note that these are 
particles that composed subhalos with more than 100 par- 
ticles that were accreted since z « 2: they make up roughly 
20% of the halo's z = particles (the other 80% were either 
added to the halo through "ambient" accretion or in subha- 
los below our particle limit or are still bound to subhalos). 
Around 75% of these particles have become unbound from 
their original subhalo due to stripping processes and are, by 
z — 0, bound only to the main halo. 

We note that unlike the background material, these par- 
ticles are not uniformly dis tributed in the halo according to 
e.g. a Navarro et all (| 19971 ) profile - or in fact according to 
any axis-symmetric profile. Instead they appear to define 
a thick plane embedded in the roughly spherical material. 
This is our third result: particles stripped from subhalos that 
are accreted from a preferred direction occupy, at z = 0, a 
confined ellipsoidal sub-volume embedded within the halo. At 
a given radial distance, the long axis of this ellipsoidal sub- 
volume is roughly aligned with that of the background dark 
halo within the same radius. 

We wish to emphasize that both the satellite infall pat- 
tern (shown in the top panels of Fig. [2| and the z = 
stripped dark matter (shown in the bottom panels of Fig. f2] 
and in Fig. [3| is remarkable constant, with respect to our co- 
ordinate system, across a large radial distance. From Fig. [3] 
we note that its is remarkable how little the orientation of 
stripped dark matter changes over such a huge distance. As 
mentioned in section 13.2.21 we wish to emphasize that the 
satellite infall pattern at 5 r v i r (roughly 1 Mpc) is coherent 
with that at distance of > 0.1 r v ; r (roughly 40 kpc). This 
brings us to our fourth main result: the angular distribution 
of satellite infall points on concentric shells centered on the 
main galaxy, is the same regardless of at what distance the 
accretion event is measured. Additionally, the dark matter 
ellipsoid composed of particles stripped from these satellites 
is oriented in the same direction across an order of magni- 
tude in extent. In other words, the 2 = position of mate- 
rial stripped from satellites that crossed 5 rvir is in roughly 
the same configuration as the 2 = material stripped from 
satellites that penetrated the inner most part of the halo, 
0.1 r v ir- Likewise, we can see a correspondence between the 
angular distribution of satellites that fall into shells of dis- 
tance > 0.2 r v i r as those that enter 5 r v i r . Note however, 
that there is less material in the inner most shells for two 
reasons: first it is a much smaller volume, but more impor- 
tantly, far fewer satellites cross the inner shells. Crossing 
the inner shells of the halo is rarer since fewer su bhalos are 
on orbits which bring them to the inner parts - 1 Gill et all 



© 0000 RAS, MNRAS 000, 000-000 



8 Libeskind et al. 



stripped DM within 0.1 r, stripped DM within 0.5r . stripped DM within 2r . stripped DM within 4 




stripped DM within 0.2r, stripped DM within 1r, stripped DM within 3r . stripped DM within 5r , 




high density 



low density 



Figure 3. The projected dark matter density ellipsoid within 8 concentric shells of radius 0.1, 0.2, 0.5, 1, 2 3, 4 and 5 r v ; r , rotated into 
the principal axis calculated at 1 r v i r . Note how the orientation of the dark matter ellipsoid's over 5 r vl[ changes little. The dark matter 
stripped from satellites accreted across 0.1 r v ; r has an orientation consistent with that material accreted at a distance of 5 r v ; r : a factor 



of 50 times greater. Note that the (bottom) panels displaying stripped dark matter within 0.2, 1, 3, and 5r v 
the same material shown in the bottom panels of Fig. 2. 



arc simply rotations of the 



found that the pericentre distribution of subhaloes 
peaks at about ~ 0.35 r v i r . Furthermore, even if a subhalo 
does penetrate these depths, it has most likely lost enough 
particles (due to tidal stripping) that it fails to make our 
mass cut. 



The degree of flattening of the two triaxial sub- volumes 
can be seen in Fig. U where we plot the axis ratio of the 
ellipsoid as a function of radius. Note that the axes ratios 
for both halos at r = r v i r indicates that these particles have 
a triaxial (i.e. a > b > c) distribution. Both halos show that 
the particles stripped from preferentially accreted subhalos 
inhabit prolate forms within r v j r , becoming more oblate out 
to r v i r - Note however, that in Fig. the reference frame is 
fixed at r v j r , indicating that even thought the ellipsoid is 
stretched along the main axis from oblate to prolate, the 
principal direction does not change, regardless of where it is 
found. 



3.3 Quantifying the preferential infall of satellites 
and their stripped dark matter 

In this section we wish to quantify what is visually apparent 
in Fig. [2] the non uniformity of satellite infall points, the 
anisotropic distribution of z = dark matter stripped from 
preferentially accreted satellites and the self-similar angu- 
lar distribution of these two components across an order of 
magnitude in radial extent. 



3. 3. 1 The clustering of satellite infall points and 2 = 
dark matter particles 

We begin to quantify this non-uniformity by computing the 
(auto) angular-correlation function of the two distributions 
shown in Fig. [21 For each shell at a given radius, we know the 
distribution of data points: the satellite infall positions and 
2 = dark matter particles. We then construct a random 
sample with the same number of points as the data. For each 
data point, we calculate the distribution of angles (9) to all 
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Figure 5. The angular auto-angular correlation function for satellite infall points (above) and dark matter stripped from accreted 
satellites (below) for 4 shells of radius 0.2, 1. 3, and 5 r v ; r (left to right) for the MW (black) and M31 (red). The blue line indicates a 
random distribution on the sky. 




M31 b/o 
MW b/a 
M31 c/o 
MW c/o 



Figure 4. The sphericity (as measured by the axes ratios) of the 
stripped dark matter particle ellipsoid as a function of radius. We 
show this quantity for the MW halo (red) and for the M31 halo 
(black). We show the the ratio of the intermediate to long axis 
(b/a) as solid lines and the ratio of the short to long (c/a) in as 
dotted lines. 



other data points DD{9), and to all random points DR(6). 
We then construct the auto-correlation function: 



w(6) = 



DD{6) 
DR{8) 



(2) 



In Fig. [S] we plot the angular correlation function of 
satellite infall points (top panels) for the same four concen- 
tric shells shown in Fig. [2] namely, r c = 0.2, 1, 3 and 5 r v i r 
for the MW (black) and M31 halos (red). Note that a ran- 
dom angular distribution would have DD(9) — DR(6) and 
thus be represented in Fig.[5]as a horizontal line at null. The 
top panels of this figure quantify what is visually apparent 
from Fig. [2] satellites that enter shells greater than the virial 
radius do so anisotropically - from preferred directions. The 
angular correlation function of satellite infall points across 
shells deep within the virial radius, is consistent with ran- 
dom. 

Like the satellite infall points, the stripped dark matter 
is far from being uniformly distributed in the halo (or on the 
sky). We show the projected angular correlation function of 
dark matter particles stripped from accreted subhalos in the 
bottom panels of Fig[3]across four concentric shells of radius 
r c = 0.2, 1, 3, and 5r vir for the MW (black) and M3 halos 
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Figure 6. The angular power spectrum of the satellite infall 
points (dot-dashed lines) with a mass weighting (dashed) and 
stripped DM (solid lines) for the MW (black) and M31 (red) ha- 
los at r c = r v i r . The plot shows the power per harmonic number, 
£, normalized by the monopole term, Co- 



(red). Within r v i r we see a clear bipolar signal indicating a 
smooth ellipsoidal distribution. As we go to larger radii, the 
bipolar signal is weakened as the dark matter becomes more 
clustered. 

We have also performed the following KS test on the 
distribution of satellite infall points. For each spherical sur- 
face we place many random strips on the sky with an angular 
thickness of 15°. For each strip we calculate the Kolmogorov- 
Smirnoff (KS) probability that the distribution contained 
within it is consistent with a random distribution. We then 
take the average KS probability of all strips to get a feeling 
for how random the sample is. We find the same conclusion: 
within r v ; r the KS test accepts the hypothesis that satel- 
lite infall points are consistent with a random distribution, 
while outside r v h-, this hypothesis is rejected at high confi- 
dence levels. Performing the same KS test on the stripped 
dark matter particles indicates an inconstancy with a ran- 
dom sample at high significance levels. 



3.3.2 Spherical Harmonic Decomposition 

The angular distribution is further investigated here by 
mean s of a spherical harmonic decomposition (e.g. IPeeblesI 
1 19801 ) of the angular distribution of satellite infall points and 
stripped dark matter particles. Although technically speak- 
ing only the satellite infall points sit on the surface of a 
sphere, in order to properly deal with the dark matter par- 
ticle distributions, we project their positions onto the surface 
across which their host satellites were accreted. 

Consider the angular distribution of N discrete objects. 
The spherical harmonic transform of the angular distribu- 
tion of these points is then given by: 



a? = jjY^W^) ( 3 ) 

i=l 

The power spectrum as a function of spherical har- 
monic, I can be calculated as 

<*=2£TlX> m | a W 

in — — £ 

The angular power spectrum is important because it tells us 
which £ mode contains the most power - or which I mode is 
contributing the most to the spherical harmonic decompo- 
sition. 

Given a dark matter particle or satellite infall point with 
coordinates (0, <j>) we may evaluate the spherical harmonic 
function for a given I and m. The original discrete point dis- 
tribution can be approximated by a smooth angular surface 
density constructed by summing over all m's and I's up to 
some ^ max according to the following equation: 

fmax 1 

a{e,4>) = Y J < Y ^( 6 ^) ( 5 ) 

f=0 m=-l 

Some of the modes of the spherical harmonic decompo- 
sition are particularly interesting. While the I — monopole 
term is trivial (in our formulation it is a constant Co = 
l/47r), the £ = 1 dipole mode and the I — 2 quadrupole pro- 
vide us with preferred directions. The relative importance of 
each mode can be obtained by computing the power spec- 
trum, Ci, according to Eq.0] 

In Fig. [II] we show the angular power spectra for the 
r v ir shell of the MW and M31 (black and red lines respec- 
tively) calculated using equation 0] Note that we normalize 
our angular power spectra by the trivial monopole (I = 0) 
contribution such that the power at a given I is measured 
relative to the monopole term. The satellite infall angular 
power spectra are shown as dot-dashed lines while the an- 
gular power spectra for dark matter particles stripped from 
preferentially accreted satellites are shown as solid lines. We 
highlight that in all four angular power spectra the maxi- 
mum occurs at I = 2, indicating that the quadrupole mo- 
ment dominates the angular power spectrum. 

In the case of satellite infall points, we note that al- 
though the I = 2 quadrupole moment is the maximum of the 
power spectrum, other moments are not necessarily negligi- 
ble. In order to further investigate the importance of other 
moments, we weigh each spherical harmonic transform by 
the mass of the accreted satellite, rrii and divide by the to- 
tal mass accreted in all satellites, M 

1 N 

a T = J-jr^ Y m i Y e m (8i, <t>i) (6) 

The mass weighted power spectrum is shown in Fig. |S]as the 
dashed lines. Weighing each infall point by the mass of the 
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Figure 7. The first three harmonics (up to the i = 2 quadrupole) of the spherical harmonic decomposition of the angular distribution 
of satellite infall points (top) and the z = dark matter stripped from these satellites (bottom) across our four fiducial shells. The blue 
dots indicate the direction of the largest eigenvector of the quadrupole tensor. 



satellite drastically increases the importance of the I = 2 
quadrupole term. This implies that although all satellites 
are being anisotropically accreted, the signal is strongest 
for the most massive satellites. This makes intuitive sense: 
the most massive satellites are born in the densest parts of 
the filament while smaller halos are more widely scattered 
across th e thick filaments, reflecting their wea ker clustering 
strength (jCole fc Kaiserlll989l : iMo et al.|[l999h . 

The dipole (and higher order moments) is relatively 
weak compared with the quadrupole. We thus focus our 
analysis on the quadrupole and its principal direction. In 
particular we can define the (dimensionless and traceless) 
quadrupole tensor, normalized by the number of particles, 
by writing 



N 



N r 2 

i=l 



Sufi ' 

~3~' 



where xf and x? are just the (x,y,z) positions of the i th 
particle. In practice the individual elements of the symmet- 
ric quadrupole tensor can b e computed directly from the 
y/ n 's (e.g. see I Jackson! 1 19621 ). Once we have obtained the 



quadrupole tensor, we can diagonalize it to obtain the princi- 
pal eigenvectors which quantify the preferred directions that 
are visually apparent in. Fig. [5] The eigenvalues are num- 
bered in decreasing ordered such that the first eigenvector 
corresponds to the largest eigenvalue (hereafter referred to 
as the "largest eigenvector"). 

In Fig. [7| we show the i — 0, 1, 2 spherical harmonic de- 
composition of the satellite infall points (top) and the dark 
matter particles contained in preferentially accreted satel- 
lites (bottom) as a function of infall radius for the MW halo, 
noting that the M31 halo looks very similar. The two direc- 
tions (positive and negative) of the largest eigenvector of 
the diagonalized quadrupole tensor (a sa tquad for the satel- 
lite infall points and cidmquad for the dark matter particles) 
are shown as the two blue dots in each panel. Note that the 
small misalignment between the two lobes of high spherical 
angular over density and the directions of the largest eigen- 
vector of the quadrupole tensor, is indicative of the existence 
of other albeit minor moments. 

When looking at the satellite infall panels, we note that 
the quadrupole dominates out to the final 5 r vlI panel where 
it begins to become weaken. Note also that although the di- 
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rection of the a sat quad changes across 5 r v i r , it does so slowly 
and incrementally. Furthermore, the location on the sky of 
fisat quad is close to that of fidm quad which also changes very 
incrementally across 5 r v i r . Note that the angular distribu- 
tion reconstructed from the first three harmonics shown in 
Fig.[7]looks qualitatively similar to Fig. [2] The largest eigen- 
vector of the quadrupole tensor points close to the "single" 
preferred direction (the red circle in Fig. [2]) identified as the 
patch on the sky with the most infall points within it. Had 
we wanted to fully reconstruct the angular surface density, 
we would have to consider harmonics beyond 1 = 2. 

We now quantify the orientation as a function of radius 
of the largest eigenvectors of the quadrupole a sat qua d and 
&dm quad in Fig. [8] for the M31 halo (top, panel a) and the 
MW halo (bottom, panel b). We show the angle between 
a S at quad and a dm quad as a function of radius as the solid 
line. For both halos, the cosine of this angle is ~ 1 within 
r v ir, indicating that the two are well aligned with each other. 
Beyond r v i r the alignment weakens for the MW at r > 3 r v h-. 
In the MW's case, the weakening of the alignment at greater 
radii is indicative of a less collimated in-flow. Stripped DM 
brought at greater radii is spread over a larger angular ex- 
tent indicating less radial satellite trajectories. For M31, the 
mild weakening of the alignment indicates that the collima- 
tion is more robust - satellites fall in on relatively radial 
trajectories. 

We now address how the largest eigenvectors of the I = 
2 quadrupole changes with respect to the virial values. The 
dotted and dashed lines in Fig. [8] show the (cosine of the) 
angle between a sat quad and a dm quad and the a sat q u ad (r vi r) 
and a dm q u ad (r v i r ), the eigenvectors at the virial radius. By 
definition both these lines lines go through the point (1,1). 
For the MW halo, (Fig. [8^) we see that these are remarkable 
stable and are well aligned out to 5 r v h-. We see a very similar 
picture for M31. 

The spherical harmonic decomposition presented in the 
preceding section allows us to quantify the visual impres- 
sion immediately apparent from Fig. [2] The principal direc- 
tions are identified as the largest eigenvectors of the 1 = 2 
quadrupole tensor. We see that indeed they point towards 
the regions of high satellite infall point density and high dark 
matter projected density. We can thus measure the orienta- 
tion between the eigenvectors and can determine how they 
move about with increasing radius. 



4 SUMMARY AND CONCLUSIONS 

In this paper we have used constrained cosmological simu- 
lations of the formation of our Local Group to examine the 
accretion of satellites galaxies. Our simulations reproduce 
the essential features of the local group allowing us to draw 
conclusions that are directly applicable to the observed Lo- 
cal environment. 




1 2 3 4 5 



r/r 

vir 

Figure 8. The orientation of the first eigenvector (a) of the 1 = 2 
quadrupole for M31 (top, panel a) and the MW (bottom, panel 
b). The solid line shows the angle between the largest eigenvec- 
tor of the satellite infall point distribution a sa t and that of the 
dark matter ellipsoid function of radius. The dotted line 

shows the value of the angle between a sa t at a given radius and 
a S at(r v ir), a * the virial radius. The dashed line shows the same 
quantitiy for the Dark matter ellipsoid, namely the angle between 
a dm and a dm (r vir ). 



We have examined the locations on the sky from which 
satellite galaxies are accreted onto the simulated MW and 
M31. We use the large scale structure to fix a coordinate 
system in which we measure the distribution of infall points 
across a number of radii. We found five important results: 

• Integrated over the entire history of the halo, subha- 
los do not enter concentric spherical surfaces (including the 
virial radius) from all directions, rather they tend to be ac- 
creted from prefered directions on the sky. 

• Part of the dark matter bound to these subhalos at 
accretion gets stripped and becomes bound to the main halo. 
At 2 = 0, however these particles define an ellipsoidal sub- 
volume embedded within the halo. 

• The orientation of this sub-volume is consistent with 
the preferred direction of the accreted progenitor subhalos: 
the preferred direction is perpendicular to the short axis of 
this ellipsoid. The dark matter stripped from these subhalos 
thus retains a memory of the prefered direction of infall. 

• The infall pattern across an order of magnitude in radial 
distance indicates that there is a coherent distribution from 
roughly 1 Mpc (5 r vlI ) to 20 Kpc (0.1 r v i r ). Satellites cross 
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concentric shells at roughly the same regions on the sky. Also 
the z — stripped dark matter in these concentric shells 
displays a remarkable angular association with the locus of 
satellites infall points. 

• By performing a spherical harmonic decomposition of 
the satellite infall points and the z — dark matter stripped 
from them, we have identified the 1 = 2 quadrupole moment 
as the dominant mode. The satellite infall point quadrupole 
moment is even more dominant when the spherical harmonic 
decomposition is weighted by the mass of the satellite. 

• The first eigenvector of the quadrupole moment, cor- 
responding to the highest eigenvalue, defines a line along 
which the satellites are preferentially accreted and where at 
z — most of the stripped dark matter particles reside. For 
the satellite infall points, this line points to within ~ 30° of 
the patch on the sky with the highest density of satellite in- 
fall points. The principal axes calculated for different shells 
are remarkably well aligned over a wide radial extent. 

Our results are important for a number of studies. 
Firstly, assuming that subhalos harbor luminous galaxies, 
and that these galaxies can be associated with the tidally 
stripped dark matter subhalos, our results confirm early 
works that attempt to explain the origin of a flattened distri- 
bution of satellites. Since the stripped material of accreted 
subhalos defines a flattened distribution at z = 0, so too 
should the galaxies they once harbored. 

The current analysis is focused on a DM only con- 
strained simulation of the Local Group, focusing on the role 
of satellites and the DM carried by the satellites into the 
main halo. This DM component constitutes roughly 20% of 
the total mass of the halo. Such a relatively small fraction 
might raise doubts on the significance of our findings on 
the distribution of that component. However, current un- 
derstanding of galaxy formation points towards a picture 
in which the galactic stellar halo is pre dominantly mad e 
of stars formed in the accreted satellites (|Bell et alj|2008h . 
Given that stars behaves dynamically in a manner similar 
to simulated DM particles, it follows that our conclusions 
concerning the DM carried by satellites pertains equally to 
the stellar halo. Our results therefore predict that the tidal 
remnants of stripped satellites (for example, Globular Clus- 
ters) should be embedded in of streams dark matter and 
stars. 

We also find a greement with stud ies of larger systems 
(e.g. clusters, (see iKnebe et all |2004 )) which have shown 
that satellites are accreted by clusters in preferred directions 
that coincide with filaments. Our results show that accretion 
of subhalo from preferred directions occurs on the galactic 
halo mass scale as well. 

Additionally, our results indicate that for galactic ha- 
los, there appears to be conherence in the infall patterns. 
Although only small percentage of satellites are on purely 
radial orbits, collapsing filaments have well defined direc- 
tions, giving rise to satellites that can penetrate deep into 



the halo's inner parts from the same directions from where 
the filament originated. 

An interesting observation can be made regarding to 
the = or ientation of th e two halos with respect to the 
larger scales. iNavarro et al.l (12004 ) argued that the rotation 
axis of a central (disc) galaxy aligns itself perpendicular to 
the minor axis of the surrounding matter. This can perhaps 
be seen in the real MW, as the North Galactic Pole points 
almost directly to the Virgo cluster. In our simulations we 
find a similar set-up: the angular momentum axis of our two 
halos seem to lie close to the plane defined by our coordinate 
system and point in the direction of our virgo cluster. 

Another interesting observation is the existence of a 
new satellit e population known as "satellites of satellites". 
Specifically, iBelokurov et al.l <|2009D discovered a new satel- 
lite of t he Milky Way - Segue 2 - embedded in a stream of 
matter, feelokurov et al.l l|2009h concluded that the stream of 
matter is the remnant of a tidally stripped parent satellite to 
which Segue-2 once belonged. Since we have found an asso- 
ciation between stripped satellite material and the location 
of satellite infall points, the fact that this object is embed- 
ded in a stream of material is readily understandable, if the 
stream is composed of (dark) matter that has recently been 
unbound and is moving with the satellite. Our results indi- 
cate that many satellites should be embedded in "streams" 
defined as coherent dark matter ellipsoids. 
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